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ABSTRACT 

We present an algorithm for quickly generating multiple realizations of A^-body simulations to 
be used, for example, for cosmological parameter estimation from surveys of large-scale structure. 
Our algorithm uses a new method to resample the large-scale (Gaussian-distributed) Fourier modes 
in a periodic A^-body simulation box in a manner that properly accounts for the nonlinear mode- 
coupling between large and small scales. We find that our method for adding new large-scale mode 
realizations recovers the nonlinear power spectrum to sub-percent accuracy on scales larger than about 
half the Nyquist frequency of the simulation box. Using 20 A'^-body simulations, we obtain a power 
spectrum covariance matrix estimate that matches the estimator in Takahashi et al. ( 2009 1 (from 
5000 simulations) with < 20% errors in all matrix elements. Comparing the rates of convergence, we 
determine that our algorithm requires ~ 8 times fewer simulations to achieve a given error tolerance in 
estimates of the power spectrum covariance matrix. The degree of success of our algorithm indicates 
that we understand the main physical processes that give rise to the correlations in the matter power 
spectrum. Namely, the large-scale Fourier modes modulate both the degree of structure growth 
through the variation in the effective local matter density and also the spatial frequency of small- 
scale perturbations through large-scale displacements. We expect our algorithm to be useful for noise 
modeling when constraining cosmological parameters from weak lensing (cosmic shear) and galaxy 
surveys, rescaling summary statistics of A^-body simulations for new cosmological parameter values, 
and any applications where the influence of Fourier modes larger than the simulation size must be 
accounted for. 

Subject headings: methods: A^-body simulations - cosmology: cosmological parameters - cosmology: 
large-scale structure of the universe 



1. INTRODUCTION 

The matter power spectrum is a sensitive probe of 
cosmological models. Observations of galaxy clustering, 
weak gravitational lensing, and the Lyman-a forest are 
all tracers of the matter power spectrum and have to- 
date provided key constraints on the initial conditions 
of the universe as well as on the growth and expansion 
history. Many future surveys will rely on tracers of the 
matter power spectrum to constrain models for dark en- 
ergy, dark matter, and massive neutrinos. However, due 
to nonlinear gravitational evolution of the matter distri- 
bution, both the mean matter power spectrum as well as 
its variance and correlations between bands in wavenum- 
ber are difficult to predict without expensive numerical 
simulations. 

Considerable work has been done to predict the mean 
nonlinear matter power spec tra to the precision necessary 



for near-ter m surveys (e.g., Smith et al.||2003 
et al. 2010). In order to use observations of t 
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is necessary to know the error distribution for the power 
spectrum estimates in addition to the mean spectrum 
predicted by the cosmological model. It is by now well 
established that nonlinear gravitational evolution creates 
a complicated error distribution for the power spectrum 
with signficantly increased small-scale variance and cor- 
relations between a ll band powers when compared with 
the Gaussia n case (|Scoccimarro et al ."1999: Meiksin &j 
' White,,1999l|Hamilton et al.l li^OO^II Neyrinck ct aL 200F 



Neyrinck fc Szapudi||2007j jfakahas hi et al.||2D09f .'"~Tlie 

ditterences between the Gaussian error and nonlinear 
error models for the matter power spectrum are quite 
significant for parameter estimation from galaxy sur- 
veys (Hamilto n et al., 2006j Neyri nck fc Szapudi 2007) , 
but uncertainty in the nonlinear galaxy bias with respect 
to the dark matter may dominate the inferred parame- 
ter constraints in the near future. Weak lensing surveys 
probe the matter distribution directly and are thus in- 
sensitive to biasing. The line-of-sight projection in weak 
lensing reduces the nonlinear contribution to the power 
spectrum covariance, but it is sti ll a significant consid 
eration for parameter estimation (Semboloni et al. 112007 
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Takada fc Jain|[2009l [Sato et al.|[2009 |. 

Standard estimators ot the matter power spectrum co- 
variance matrix in a given model require a large ensem- 
ble of statistically independent realizations of the matter 
density over volumes at least as l arge as the survey be- 
ing analyzed. Mciksin & White p999j found that at 



least several hundred realizations are necessary to esti- 
mat e the covariance over a n interesting dynamical range 
and Takahashi et al. ( 2009 1 have used 5000 simula tions to 
obtain a low-nois e estimator for a Gpc"^ volume. 



Hamil- 



ton et al. (2006) attempted to estimate the covariance 
matrix by sub-sampling a single simulation volume, but 
found that the window functions they applied to define 
different pseudo- independent subvolumes significantly al- 
tered the covariance. That is, the Fourier transform of 
a windowed density field has Fourier modes that are lin- 
ear combinations of the modes in the iV-body simulation 
volume with periodic boundary conditions. The modes 
of the windowed density then have different covariance 
than the modes in the periodic box, which is particu- 
larly enhanced by t he correlations between the largest 
and smallest scales f Rimes fc Hamil ton 2006). In a sim- 

showed that 



(j2009) 



ilar line of inquiry, [Norberg et al. 

jacknife estimates oi the covariance matrix from simu- 
lated survey regions systematically underpredict the true 
nonlinear covariance matrix. Finally, if some information 
is known about the structure of the power spectrum co- 
variance, th en the "shrinkage estimator" of |Pope fc Sza-| 
pudi ( 2008 ) gives a way to reduce both the bias and noise 



m covariance estimates from a fixed number of density 
realizations. However, in ou r experience the "shrin kage 
estimators" as formulated in Pope & Szapudi (|2008 1 rely 
quite heavily on the prior information one is able to sup- 
ply- 

In this paper we describe a method to obtain multiple 
realizations of the matter density from a single iV-body 
simulation by resampling those Fourier modes (in a pe- 
riodic box) that can be approximated as Gaussian dis- 
tributed with zero correlations between Fourier modes. 
While we anticipate this algorithm could be useful for a 
number of applications (suc h as adding power larger than 

Colell997), in this paper we 
if 



the simulation volume as in 
focus on the minimum num 



ber ot simulations needed to 
estimate the matter power spectrum covariance matrix. 
As a benchmark fo r our study we use the results of Taka- 
hashi et al. (20091, who used 5000 iV-body simulations 



to estimate the covariance matrix in a 1 {h~^Gpc)^ sim- 
ulation volume. While we focus here on the 3-D matter 
power spectrum, one final goal will be to efficiently es- 
timate the weak lensing power spectrum covariance for 
upcoming cosmic shear surveys. An exa mple of t his us - 
ing a brute- force approach is given in Sato et al. ( 2009 ) . 

The structure of this paper is as follows, in Section 2] 
we describe our method for resampling large-scale modes 
in a periodic A'^-body simulation and the necessary ad- 
justments to the remaining small-scale modes. We give 
an approximate algorithm requiring only a single snap- 
shot, but limited t o simu lations of at least several Gpc on 
a side, in Section |2.1.2| We describe our full algorithm 
for arbitrary simulation sizes (as long as t he lar ge-scale 
modes are Gaussian distributed) in Section 2.1.3 In Sec- 
tion [3] we present validation tests for adding new large- 
scale mode realizations. We describe how our resampling 



algorithm can be useful for p ower spectrum covariance 
matrix estimation in Section 13.2.31 We discuss several 
future directions and applications for this work in Sec- 
tion |4j Where numerical results are given, u nless other- 
wise s t ated, we assume the same cosmology as 'Tak ahashi] 
eraLl ( [2009| with = 0.238, fli, = 0.041, Qjy = 0.762, 



0.958, as = 0.76, and Ha = 73.2 kms^i Mpc"'. 

2. RESAMPLING LARGE-SCALE POWER 

Our goal is to isolate and resample those Fourier modes 
in a single periodic A^-body simulation that can be ac- 
curately approximated as Gaussian distributed. We as- 
sume that the Fourier modes of the initial conditions for 
the simulation are generated from a random realization 
of a complex Gaussian field on a grid from a specified 
power spectrum (i.e. the real and imaginary parts of the 
Fourier modes are each drawn from independent Gaus- 
sian distributions such that the amplitude of the modes is 
equal to the square root of the power spectrum) . As the 
simulation is evolved, gravitational interactions between 
the particles tracing the matter density field cause the 
phases of the Fourier modes to become correlated and 
the distribution of Fourier mode amplitudes to become 
skewed towards larger values. The growth rate of cor- 
relations and departure from Gaussianity are functions 
of wavenumber as gravitational collapse proceeds most 
rapidly on small-scales. 

So, we first need to separate the Fourier modes into 
"large-scales" defined to be Gaussian distributed, and 
"small-scales" that we will define to be all the remaining 
Fourier modes in the simulation. Because it worked for 
our purposes we crudely separate large and small scale 
modes by defining all modes with wavenumber amplitude 
less than some value fcthrcsh as Gaussian distributed. In 
practice we defined fcthrcsh so that the r.m.s. mass over- 
density smoothed in top-ha t spheres of radius 27r/A;thresh 
was less than 0.2 (similar to Cole|1997 |. This is based on 
the idea that larger mass overdensity fluctuations indi- 
cate more nonlinear growth leading to larger departures 
from Gaussianity. 

Having defined those Fourier modes that we will model 
as Gaussian, we next resample them by generating new 
realizations of the initial^conditions that have nonzero 
Fourier modes only for < fcthrcsh and extrapolating 
the amplitudes of the modes to the desired redshift using 
linear theory. 

To obtain a new realization of the density field that 
approximates what we would have obtained had we run 
a simulation with the new large-scale modes in the ini- 
tial conditions of the simulation, we have to adjust the 
small-scale modes in some way to account for the phase 
correlations that would have emerged between the large 
and small scales. 

2.1. Resampling algorithm 

In this Section we imagine that we know the par- 
ticle positions in a simulation that has no large-scale 
power (where "large-scale" denotes Fourier modes with 

|fc| ^ fcthresh) and we seek an algorithm to add new 
large-scale modes. We will denote the simulation with- 
out large-scale power as a "small-modes" simulation and 
the simulation with all Fourier modes nonzero as an "all- 
modes" simulation. 
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2.1.1. Physical picture 

Altering the amplitudes and phases of the large-scale 
modes of the mass density in a cosmological volume has 
two dominant effects, 

1. sub- volumes are expanded or compressed, 

2. sub- volumes gravitationally evolve with an effec- 
tive matter background density that is enhanced 
or diminished, 

depending on where the large-s cale modes become more 
or less ovcrdense, respe ctively ( Frenk et al."1988( iLittle 



et al.|19 91; Tormcn & Bertsching er 1996; Cole 1997fp^ 
gulo &: White||2010| ). We follow Cole (1997) to moHeT 
both of these effects. To model the expansion and com- 
pression of sub- volumes, we perturb the particle posit ions 
using the Zel'dovich displacements ( Zerdovich|1970 1 cal- 
culated from the large-scale modes, d^. 



dL(x, a) 



(27r)3 



fc2 



(1) 



The expansion and compression of sub-volumes is anal- 
ogous to the frequency modulation of a radio wave. The 
large-scale modes to be added are a new signal that mod- 
ulates the sizes of sub- volumes in the matter density field; 
expressed to first order by the Zel'dovich displacements 
in equation ([T]). 

Compressea (expanded) regions also evolve with an in- 
creased (decreased) effective matter density, altering the 
growth rate of structure. Because the growth rate in- 
creases monotonically with time, we model the changing 
effective local matter density by finding new times a'(x) 
at which the linear growth rate matches the growth rate 
in a universe with matter density r2„i (1 + <5l(x, a)) (with 
r^A kept fixed). The perturbed scale factors a'(x) are 
then defined by the relation. 



D{a', n^) « D{a, (1 + ^l(x, a))). 



(2) 



where D{a, flm) is the linear growth function. In prac- 
tice, we solve equation ^ numerically to get a' as 
demonstrated in Figure [T] 

2.1.2. Approximate method for a single snapshot 

Combining the two effects of expansion/compression 
and time-perturbations, the model for adding large-scale 
fluctuations to the mass density field is. 



J(x, a) = (5l(x, a) + (5s(x', a'), 



(3) 



where (5(x, a) = p(x, a) j 1 is the mass over-density that 
would be measured in the all- modes simulation, (55(x, a) 
is the mass over-density measured in the small-modes 
simulation, Si^ has zero power for |fc| > fctiu-eshi is de- 
fined by equation (l2|, and x' = x -I- d2,(x, a). Note that 
for this definition of x' the Zel'dovich displacements are 
added to the particle positions prior to application of the 
time-perturbation (i.e., at the positions at time a, not at 
a'). We will discuss this choice below. 

/ 

When \a! - a| < 1 (which requires ^(^^(x)! \ < 1) 
the first order correction to 8s in a! should be sufficient 
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Figure 1. Top left: Linear growth function normalized to one 
at the time of the simulation initial conditions (aic = 1/21) ver- 
sus scale factor. Top right: Linear growth function versus matter 
density, also normalized to unity at the time of the initial condi- 
tions for all values of Q,m- Bottom: scale factor as a function of 
matter density at matched linear growth. The vertical dashed line 
shows the value of = fi^'^j = 0.238 used in our simulations. The 
vertical dotted lines show the values of it 4.3a-(<5i)), where 

""(^l) ~ 0.17 is the r.m.s. mass over-density in top-hat spheres at 
the scale R = 2tt / {25fik p) (kp = 27r/Li,ox) for the density per- 
turbations with nonzero Fourier amplitudes only for |fc| < ^thresh- 
The horizontal dotted lines in the bottom panel show the corre- 
sponding scale factor range over which our small-modes simulation 
must be saved. 



for adding 5tAo the density field. Expanding both sides 
of equation (pi) to first order in (a' — a) and Sl, 



D{a') w D{a) + (a' - a) 



dD 
da 



(4) 



D{a, r2,„(l + (5L(x,a))) « D{a, $!,„) -I- fi„(5L(x, a) 



dD 



Solving for a' gives, 

1 « (5L(x,a) 



a 
a 



d\nD d\nD 
dlnflm din a 



(5) 



which shows that the perturbed scale factor is roughly 
proportional to the large-scale density perturbations be- 
ing added. Using equation ([s]) in the Taylor expansion 
of equation ^ gives, 

(5(x, a) w Sl{x, a) + (5s(x', a)+ 

.d6s{^',a) fdlnD d\nD\ 
^L[^,a)^r. hjl— Ft/t71 ■ (6) 



9 In a 
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This suggests that the addition of large-scale modes can 
be quickly computed given an estimate of the nonlinear 
growth rate of the small-modes simulation. With a sin- 
gle snapshot, the growth rate can be estimated from the 
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continuity equatioij^ 

(9(5(x, a) 



H{a)- 



9 In a ( 
and, using equation ([s]), 

d5s{'x.,a) 9(5(x, a) 



-V.((l + J)v) 



9 In a 



9 In a 



din a 



Jl(x, a) 



(7) 



(8) 



So, given particle positions and velocities at scale factor 
a for a simulation that has no large-scale modes, equa- 
tion ^ provides a method to add new large-scale mode 
realizations when the large-scale density perturbations 
are small. 

According to equation ([5]) the range of a' — a is pro- 
portional to the range of Ti(x,a). Because (5l(x, a) is 
Gaussian distributed by construction, we can determine 
the range of 5^ from its variance. 



((5L(x,a)) 



d\nkA\k)Wl^_^^,{kR), (9) 



where kp = 27r/Lbox is the fundamental frequency of the 
simulation box, A^(fc) = k^ P{k) / {2'it'^) is the theoretical 
dimensionless matter power spectrum, and Wtop-hat(fc^) 
is the Fourier transform of a spherical top-hat window of 
radius R. The scale R is roughly given by either the 
pixel scale of the gridded density field or the mean inter- 
particle spacing, whichever is larger. For the simulations 
we use below, Lbox = 1000 /i~^Mpc, fcthrcsh — ^kp w 
0.05/i~iMpc, and, with 256^ particles, R - 2Ti/{2^Qkp). 
Evaluating equation ^ gives a « 0.17 which is still less 
than 0.2 as we required when choosing fcthresh- But, as 
shown in Figurejl] the variations in 5l(x, a) greater than 
about 2(7 lead to|a' — a| > 1 making equation Q insuf- 
ficiently accurate. In our test simulations Sp actually 
varies over ±4.3(t « ±0.87, giving a range of perturbed 
scale factor 0.26 ^ a' < 4.4 (shown by the dotted lines 
in Figure [T]). 

In Figure [2] we show the maximum value of \a' — a\ 
as a function of the simulation box size while keeping 
kthiesh/ kp fixed at 8 and, as in the preceding para- 
graph, assuming maximal variations in 5p{'x.,a) of 4.3cr. 
For Lbox ~ 3000/i~^Mpc, as in the recently completed 
MXXL simulatiorj^ the variation in a' — a is less than 
0.2 and the approximate formula in equation ([6]) could be 
sufficient for adding the large-scale modes to the small- 
modes simulation. To further assess the validity of equa- 
tion ([g]), we show a histogram of the values of the third 
term in equation ([6|, normalized by (5s(x', a), on a 512'^ 
grid in our 1 /i~^Gpc small-modes simulation in Figurejsj 
The distribution of grid cell values is strongly peaked 
around zero, but has wide tails extending to values of 
several hundred. This again illustrates that equation (|6| 
is a bad approximation for a (1 /i~^Gpc)'^ volume. How- 
ever, if we scale the bin values on the lower horizontal 

^ Note that for simulation volumes 2> 1 Gpc on a side rela- 
tivistic corrections to the continuity equation could be important. 
These arise from the way that the large scale potential perturba- 
tions modify the geodesic equations. 

^ The MXXL simulation (Angulo et al. in preparation) was run 
by the Virgo Consortium (http://www.virgo.dur.ac.uk) in the 
spring of 2011. It is 3 h~^G\>c on a side with 3 X lO^'^ particles 
with mass 6.3 X 10^ Mq. 
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Figure 2. Maximum variation in perturbed scale factor a' as a 
function of the simulation box size for fixed fcthresh/^F = 8 and 
assuming the maximum fiuctuations in 5j^(x, a) are 4.3cr. Top: 
Maximum absolute difference between a' and a = 1 for the growth 
matching condition for adding large-scale modes to a simulation 
(equation [2]l . Bottom: Maximum a' for the growth matching con- 
ditio n f or removing large-scale modes from a simulation (equa- 
tion |11^ . There is no solution for a' left of the vertical dashed 
line. 
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Figure 3. Histogram of the values of 



{a' — a) din 



.512^ ! 



— n til J. ^ ijnd 

a In a 

for our ACDM cosmology in a 1 (h^^Gpc)"^ simulation. The tick 
marks on the top (blue) axis show the projected bin values if we 
scale the bins on the bottom axis by the ratio of the maximum 
perturbed scale factor, max(a' — a), in a 3 h~^Gpc box versus 
the maximum perturbed scale factor in our 1 /i-^Gpc box. The 
vertical dashed lines show the values of ±1 on this axis, where the 
third term in equation (jcjl would be small relative to Sg. 
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axis by the ratio of the maximum perturbed scale factor 
max(a' — a) in a (3 h~^Gpc)^ simulation to the maximum 
(a' — a)in our simulation, we get the axis on the top of 
Figure Is] For a (3 h~^Gpc)^ simulation 92% of the grid 
cells oi^he third term in equation ^ have values less 
than one. So, equation ^ may be a good enough ap- 
proximation for adding large-scale modes to the MXXL 
or other simulations with similarly large volumes. 

2.1.3. Full recipe 

For box sizes smaller than several Gpc on a side, 
we need an algorithm that accurately implements the 
growth matching condition of equation (pi) to compute 
the density 6s{x',a') in equation ([3]). F'or co mputa- 
tional expediency and for later comparison with |Taka- 



hashi et al. 



( 2009 ) , we have limited our test simulations 
iOOO/i~^Mpc. Therefore, we have not tested 
Instead, because our test 
> 1, we apply the 



to Lbox 

the accuracy of equation (6l 
simulation parameters yield 

model of equation ([3]) by adjusting the positions of each 
particle in the simulation. If Xi(a) is the position of the 
ith particle at scale-factor a in the small-modes simula- 
tion, the model from equation ([3| with modes (5^ added 
to the simulation can be rewritten for each particle po- 
sition as. 



X- {a';SL) = [xj(a) -l-di (x,(a), a*)]^^^, , 



(10) 



where a* is the scale factor of the target time where 
the full density field (with new large-scale modes) is 
desired. That is, we first apply the same Zel'dovich 
displacement to the particle positions in all the snap- 
shots (evaluated at the unperturbed positions in each 
snapshot). Then, we interpolate between these modi- 
fied snapshots to find the particle position at perturbed 
time = a'{xi + d((5i(xi(a)), a*)). To find the parti- 
cle positions at all values of a' we must save simulation 
snapshots at times covering the range of a' values and 
spaced sufficiently closely in a so that the interpolation 
of positions is accurate. This requires running the simu- 
lation far into the future (a > 4 for our test simulations) 
and saving ^ lO's of snapshots, which creates both an 
extra computational and also a disk storage burden over 
the approximation in equation Q. 

The steps of our full algorithm for adding large-scale 
modes to the small-scale density are then: 

1. Choose a threshold scale fcthrosh defining the bound- 
ary between large and small scale modes and a tar- 
get time a* where the resampled density is desired. 

2. Calculate the range of perturbed scale factor a' de- 
fined by equation ([2]) about a* where the simulation 
snapshots will be required. 

3. Run an iV-body simulation in a periodic volume 
with nonzero Fourier modes in the initial condi- 
tions only for |fc| > fcthrcsh and save snapshots over 
the required range of a' that are spaced sufficiently 
close together to allow for accurate interpolation of 
particle positions. 

4. Calculate a new realization of large-scale modes 
i5i(a*) on a grid. 



5. Calculate the Zel'dovich displacements defined in 
equation ([T]) from S^. 

6. Apply the Zel'dovich displacements to each snap- 
shot of the small-modes A^-body simulation. 

7. Move each particle in the small=modes simulation 
to its location at time a'(x(a) + d^) as defined by 
equation [2] interpolating between snapshots when 
necessary. 

The resulting particle positions give a close approxima- 
tion to the particle positions that would be obtained by 
running an A-body simulation with the Fourier modes in 
Sl present in the initial conditions (with the amplitudes 
suitably scaled according to the linear growth function). 

The two operations of adding Zel'dovich displacements 
and perturbing the time where the particle positions are 
evaluated do not commute when the density field has 
evolved nonlinearly from the initial conditions. To de- 
cide on how to order these operations we imagine an ap- 
proximation to an A^-body simulation where the initial 
conditions are first evolved by adding Zel'dovich displace- 
ments to the particles and then the density is scaled by 
the linear growth function. Reversing these operations, 
with the Zel'dovich move applied last, would be less anal- 
ogous to the way the particle positions actually evolve in 
a simulation. In this latter scenario, the particles would 
first gravitationally collapse at rates determined by the 
local mean matter density and then at late times would 
undergo displacements to expand or compress regions in 
over- or overdense volumes. There is nothing physical 
about this ordering of operations. We have confirmed 
with our test simulations that the mode resampling does 
not work well unless the Zel'dovich move is the first step 
of the mode-resampling algorithm. 

2.2. Algorithm implementation 



The algorithm in Section |2.1| describes how to add 
large-scale Fourier modes to a simulation that had zero 
large-scale power in the initial conditions. Ideally we 
would like to first remove the existing large scale modes 
from a previously run simulation and then apply the al- 
gorithm in Sectio n [2] We will discuss why this is chal- 
lenging in Section 2.3 but first we describe details of the 
implementation of the mode addition algorithm. 

To generate a simulation without large-scale power 
(what we call a "small-modes" simulation), we gener- 
ated Gaussian Fourier mode realizations on a grid only 

for those modes, set all modes with |fc| < A:thrGsh to zero, 
and then performed the reverse Fourier transform to get 
the initial density field that is used for the first particle 
displacement step. We then evolve the simulation nor- 
mally, but saving 43 snapshots equally spaced in ln(a) 
over the range 0.1 < a < 5.4 so the particle positions 
at perturbed times can be accurately determined by in- 
terpolating between snapshots. For testing purposes we 
ran another simulation using the same initial conditions, 
but before the large-scale modes were set to zero (called 
an "all- modes" simulation). The particle positions at 
a = 1 in the all-modes simulation give our benchmark for 
testing the reconstruction of the density field when the 
large-scale modes that were deleted from the initial con- 
ditions are added back into the small-modes simulation 
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n,„ =0.133 
n,„ =0.306 
Full box 
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Figure 4. Matter power spectra for several snapshots of the sim- 
ulations with scale factor increasing from bottom to top. The 
solid lines are for the all-modes simulation, dashed lines are for 
the small-modes simulation and dotted lines show the linear the- 
ory prediction. Notice the divergence of the nonlinear all-modes 
and small-modes spectra as the simulation proceeds. 

at a = 1. Note that all the phases in the initial conditions 
were identical for both the all-modes and small-modes 
simulations. Unless otherwise specified, all our simula- 
tions were in box 1 h~^Gpc on a side with 256"^ particles 
and were evolved using a lean version o f the Gadget-2 
code ( jSpringel et al.|2001 Springel|2005 ) from Gaussian 
initial conditions at redshiit 21). 

Because gravitational evolution couples small-scale 
Fourier modes to the largest modes in the simulation 
box, the nonlinear perturbation growth rate in the small- 
modes simulation (i.e., with large-scale modes missing) 
is smaller than that in the all-modes simulation, as illus- 
trated with the power spectrum comparison in Figure |4j 
The solid lines show the power spectra of select snapshots 
in the all-modes simulation while the dashed lines show 
power spectra of the small-modes simulation at the same 
times (recall that the power spectra are equal in the ini- 
tial conditions). At late times, the all-modes power spec- 
trum has a significantly larger amplitude on small scales 
than the small-modes power spectrum. A first test of 
our mode-addition algorithm will be to correct for this 
deficit of nonlinear power. 

In order for the growth-matching condition in equa- 
tion ([2]) to provide a useful means for correcting for 
missing large-scale power, the errors in matching the 
growth rate by perturbing the matter density must be 
much smaller than the change in the growth rate due 
to the ommission of large-scale power. We tested this 
condition explicitly by running two simulations with dif- 
ferent values of and we plot their power spectra at 
matched linear growth in Figure [5] The power spectra 
in the two cosmologies match to < 0.1% for wavenum- 
bers k < 0.2 /i/Mpc, which is much less variation than 
between the all-modes and small-modes power spectra. 



Figure 5. Power spectra relative to the all-modes spectra from 
simulations with different values of Qm plotted at matched values 
of the linear growth. For comparison, the power spectra for our 
fiducial cosmology with and without large-scale modes are also 
shown. The variation in the spectra with different Qm is only 
apparent for k > 0.5 hMpc~^ and even there it is much less than 
the difference between the spectra with large-scale modes removed. 

This is why matching the linear growth is a good way to 
mim ic the change in th e local matter density. We note 
that Zheng et al. ( 2002 ) have performed similar tests and 
showed that when the linear theory power spectrum is 
fixed (as we have assumed), several clustering statistics 
besides just the nonlinear power spectrum can be easily 
related between models with different f2m and cts- 

Next, we tested the vario us s teps of the mode addi- 
tion algorithm from Section |2.1| by matching the parti- 
cle positions between the all-modes and the small-modes 
simulations by adding the large-scale modes from the all- 
modes initial conditions back into the small-modes simu- 
lation at a = 1 (after scaling the large-scale modes by the 
linear growth function). Figure [6] illustrates our compar- 
ison for a section of the l(/i^^Gpc)'^ simulation volume. 
The black points show the particle positions in the all- 
modes simulation while red shows the positions in the 
small-modes simulation. It is readily apparent to the eye 
that both the time-perturbation and Zel'dovich moves 
are necessary to obtain a good match to the particle po- 
sitions in the all-modes simulation. The Zel'dovich move 
alone does a better job matching the particle positions 
in the largest over-densities than the time-perturbation 
alone, but leaves significant position offsets compared to 
the final algorithm combining both perturbations. 

2.3. Subtracting large-scale power 

The algorithm we have demonstrated requires that an 
A'^-body simulation be run with initial conditions that 
have the large-scale Fourier mode amplitudes set to zero 
and with a relatively large number of late time snapshots 
saved that can be used for interpolating particle posi- 
tions. Ideally, we would like to have an algorithm that 
can be applied to existing A^-body simulations. Such an 
algorithm would first require a method for subtracting 
the large-scale modes already present in the simulation. 
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Figure 6. Particle positions in a section of our 1 Gpc^/h"^ simulation box (of width 50 Mpc/h). The red points denote the particle 
positions in a simulation with the large-scale Fourier amplitudes omitted from the initial conditions as described in the text. The black 
points denote positions in a simulation with all Fourier modes included, but that is otherwise identical to the simulation denoted by the 
red points. The 4 panels show various components of our algorithm for replacing the missing large-scale power. The top left panel has no 
correction applied to the red points, the top right has the red particles shifted to their positions at different times so that the linear growth 
rate mimics the perturbation to the local matter density from the added large-scale modes. The bottom left panel shows a Zel'dovich move 
applied to the red points, where the Zel'dovich displacements are derived from the large-scale modes only. Finally the bottom right panel 
shows our full algorithm combining the time-perturbation and Zel'dovich displacements. The axis labels are in /i~^Mpc. 

Our procedure for adding large-scale modes should work 
in reverse, but with the new growth-matching condition, 



D{a\ + <5l(x, a)) ^ D{a, n„,). 



(11) 



We arrive at this condition by modeling regions where 
6l is overdense (overdense) as regions that have evolved 
without large-scale power but with a larger (smaller) 
local value of flm- We then find the time a' where a 
simulation with density ilm(l + <5l(x, a)) and no large- 
scale power would have the same linear growth as the 
original simulation with density flm at time a. How- 
ever, removing the effect of large-scale under-densities 
can require evolving those overdense regions far into the 
future to match the growth on the right-hand side of 
equation (111. For ACDM, the local linear growth in 



large under-densities can freeze an d will never match 
the right-hand side in equation (11) even if evolved in- 
finitely into the future. The growth matching condition 



in equation (11) and the problem with under-densities 
is illustrated m Figure [Tj For our fiducial cosmology, 
-^box ~ 1 h^^Gpc, and fcthrcsh — 8fc_F all regions in 5^ that 
have density smaller than the mean density minus ~ 2a 
will be impossible to correct using the growth match- 
ing condition of equation (11). In the bottom panel of 
Figur e ^ we show the maximum a' derived from equa- 
tion ( |11| as a function of simulation box size assum- 
ing kthrosh/kp = 8 and that the maximum variations in 
^l(x, a) are 4.3a. For cubic simulation boxes with sides 
less than ~ 1.8 h~^Gpc there is no solution for a' . For 
larger box sizes the maximum a' quickly decreases from 
~ 6 to very close to one. So for simulation box sizes of 
several h~^Gpc, removing large-scale modes may be fea- 
sible with minimal extra computational effort. However, 
again for computational expediency, we have limited our 
considerations in this paper to our 1 h^^Gpc simulations 
and therefore do not consider the removal of large-scale 
modes further. 
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Figure 7. Scale factor versus matter density at matciied linear 
growth for mode subtraction. Left panel: Tlie scale factor versus 
the linear growth function normalized to one at the time of the 
initial conditions (a/c = 1/21) for different values of Qrn over the 
range [0.01,0.5]. Each black line represents a different value of f!,„ 
starting with Qrn = 0.01 for the left-most line and increasing in 
steps of 0.05. The solid vertical red line shows the value of the 
normalized growth function at a = 1 and the target mass density 
= 0.238. Right panel: The growt h ma tching condition for re- 
moving large-scale modes in equation amounts to finding the 
value of the scale factor in the left panel where the red vertical 
line intersects the black lines for each Qm value. The line of inter- 
section points is plotted as the black line in the right panel. The 
vertical dotted lines show the range of 0^(1 it 4.3(t((5l)) (as in 
Figure The scale factor value matching the upper bound in the 
perturbation to Qm is shown by the dotted horizontal line. There 
is no scale factor value that can match the linear growth at the 
lower bound of the perturbed Qm, rendering the mode subtraction 
impossible for the overdense regions of 5^ ■ 

3. VALIDATION TESTS 

We present further tests in this section to assess the ac- 
curacy of our method for introducing large-scale modes 
in A^-body simulations. We are primarily interested in 
the matter power spectrum and power spectrum covari- 
ance matrix and investigate estimators for these in turn. 

3.1. Recovering the nonlinear power spectrum 

In Figure |8] we demonstrate the reconstruction of the 
nonlinear power spectrum from the small-modes simu- 
lation. Using the all-modes simulation as our bench- 
mark, we used the Fourier modes from the all-modes 
initial conditions with |fc| < fcthrcsh to determine the 
Zel'dovich moves and perturbed times for adjusting the 
small-modes simulation (so if our algorithm worked per- 
fectly the power spectra should be perfectly matched). 
In addition to the ACDM simulations described in Sec- 
tion |2.2[ we ran analogous all-modes and small-modes 
simulations for an SCDM cosmology (i.e., with fl„i — 1). 
Because the SCDM model has no dark energy, the linear 
growth rate is completely determined by the matter den- 
sity, and our growth matching condition in equation ([2| 
should be an especially good approximation for the effect 
of large-scale Fourier modes on the gravitational growth 
rate. Figure [8] shows that the fractional errors in the 
reconstructed power spectra are indeed slightly smaller 



for the SCDM simulation, but the errors are much less 
than 1% for both cosmologies. Due in part to the cloud- 
in-cell deconvolution (which over-corrects the power on 
small scales), the fractional errors begin to increase sig- 
nificantly at about O.SfcNyquist- 

Note that the Zel'dovich correction alone, shown by 
dotted blue lines in Figure |8j does a very poor job of 
correcting the small-scale power spectrum for the added 
large-scale modes. This seems to be in contrast to the 
lower left panel of Figure |6] where the Zel'dovich correc- 
tion alone appears to correct the particle positions rather 
well. We now see that the eye notices the agreement of 
large-scale structures more easily than the disagreement 
of small-scale structures in the dot plots of Figure |6] 

3.2. Power spectrum covariance matrix 

The gravitational growth of perturbations increases the 
small-scale power spectrum variance and induces signif- 
icant correlations b etween power s pectrum bands ( Sco(>] 
[cimarro et aH1999{|Meiksin fc White[l999 ). In thisEec^^ 
tion we test how well our mode addition algorithm can 
reproduce these effects. 

By perturbing a single small-modes simulation with in- 
dependent Gaussian realizations of the large-scale modes, 
we can generate many realizations of the density field. 
These are not independent realizations because the 
small-scale modes are initially the same for each resam- 
pling. We use the covariance matrix estimated in |Taka-| 



hashi et al. ( 2009 ) from 5000 A^-body simulations to test 
our sample covariance estimates from simulations with 
resampled large-scale modes. 

3.2.1. Sample covariance matrix estimator 

First we define the estimator for the covariance matrix 
that we use for our comparisons. 

Let 6 ft denote the amplitude at grid point fc^ of the Dis- 
cret Fourier Transform (DFT) of the density field from 
an A^-body simulation measured on a 3D uniform peri- 
odic grid. Then, the standard estimator for the power 
spectrum of the mode amplitudes is obtained by averag- 



ing 



over a spherical shell of constant radius 



(12) 



where Nh. is the number of modes that fit in the shell. 



We follow Takahashi et al. ( 2009 1 and set this shell thick- 
ness to 0.01 /iMpc~ . If we generate Nj- realizations of 
the density field (e.g. by running Nr identical iV-body 
simulations except with different random number seeds 
to generate the initial conditions), the mean power spec- 
trum estimator for all the realizations is: 



Pn 



1 



(13) 



The estimator for the covariance of the power spectrum 
from all the realizations is then: 



(14) 
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Figure 8. Power spectra relative to tiie all-modes spectra demonstrating the accuracy with which large-scale power can be added to 
a simulation without any large-scale modes. Left: SCDM, right: ACDM. The large-scale modes added to the simulations without 

large-scale power (small-modes) have the same phases as in the simulation run with all modes present (shown by the gray horizontal lines). 
The vertical dashed line marks the Nyquist frequency of the density fields. 

volution in Fouri er space imposed by th e configuration- 
space windows of Hamilton et al. (2006), we attempted 
to sub-sample the i'burier modes that enter the shell- 
averaged power spectrum estimator in equation ( 13 ). For 



iMeiksin fc White! Jl999[ ), [Norberg et "aT] ( |2009[) , and 
[Takahashi et ai. ( |2U09p toimd that at least several hun- 
ored realizations ot the matter density field are required 
to accurately estimate the sample covariance matrix. 
This means that hundreds (or thousands in some cases) 
of A^-body simulations must be run with different ran- 
dom number seeds in order to obtain a power spectrum 
covariance matrix estimate that is sufficiently accurate 
for cosmological parameter estimation. 

3.2.2. Multiple power spectra from a single simulation box 

Because running hundreds or thousands of A'^-body 
simulations is often prohibitively expensive, previous 
studies have attempted to estimate the power spectrum 
covariance from a single iV-body simulation by apply- 
ing window functions to the gridded den sity field t o sub- 
samplc the simulation volume. Hamilton et al. (20061 
found that because the windows are convolved with the 
Fourier modes of the periodic simulation box, the result- 
ing power spectrum covariance matrix estimates are sig- 
nificantly biased with respect to the result obtained from 
an ensemble of periodic simulations. Furthermore this 
bias cannot be simply accounted for by deconvolving the 
window functions from the power spectrum estimates, 
but is a result of higher order connected correl ations en- 
tering the c ovariance matrix estimator, which [Hamilton| 



et al. ( |2006[ ) call "beat-coupling." 

Although direct sub-sampling of the simulation volume 
has been shown not to reproduce the covariance from 
an ensemble of simulations, the idea that different re- 
gions of a large simulation volume provide statistically 
independent information about the small-scale covari- 
ance is sound if we accept the ergodic hypothesis and 
assume our simulation is big enough. To avoid the con- 



shells with large wavenumbers there are thousands of 
modes used to estimate the power spectrum. The non- 
linear power spectrum covariance depends on the con- 
nected four-point correlations of the Fourier modes, but 
with thousands of modes on a uniform grid there must 
be many redundant four-point configurations that could 
yield independent samples of the (shell-averaged) power 
spectrum covariance (after appropriate scaling by the 
number of modes used to estimate the mean power spec- 
trum) . 

The enumeration and counting of the four-point con- 
figurations in each shell is complicated, so we instead 
divided the surface of each wavenumbe r shell into uni- 
form area HEALPix( Gorski et al. [2005 1 pixels and then 
selected only one wavevector within each pixel to ob- 
tain a power spectrum estimate. For a given number of 
pixels covering the spherical shell, the minimum multi- 
plicity of wavevectors in a pixel determines the number 
of pseudo-independent power spectrum estimates we can 
derive with this method. That is, if all pixels have at 
least n wavevectors within their boundaries, then we can 
generate at most n pseudo-independent realizations of 
the power spectrum estimates. We defined smaller (and 
therefore more) pixels for shells of increasing wavevector 
magnitude. This ensured that the multiplicity of each 
pixel was greater than one for a wide range in wavevec- 
tor magnitudes. 

We applied this sub-sampling to a small-modes simula- 
tion with 500 large-scale mode resamplings. The result- 
ing power spectrum estimates and correlation coefficients 
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Figure 9. Top: Power spectrum estimates from using only a sub- 
sample of the Fourier modes in each wavevector shell. The power 
spectra are normalized by the power spectrum of the all-modes 
simulation at a = 1, so the wiggles at small k are indicating the 
sample variance in the all-modes simulation. Each power spectrum 
estimate from sub-sampled fe-shells is shown by the grey lines in the 
top panel. The bottom panel of the top figure shows the number 
of Fourier modes used to estimate the shell-averaged power spec- 
trum. The large jumps are where we changed the number (and 
size) of the pixels used to sub-sample each fc-shell. Bottom: Cor- 
relation coefficients of the p ower spectrum estimate s (solid, black) 
compared to the result from [Takahashi et al.| | [2009[ l (dashed, red). 

of the estimated power spectrum covariance are shown 
in Figure [9j While the mean power spectrum is suc- 
cessfully recovered with this algorithm, the power spec- 
trum covariance estimator has nearly zero off-diagonal 
terms, in stark contrast to the result from Takahashi, 



et al. (2009) (shown in red, dashed lines) for wavenum- 
bers > 0.1 hMpc~^. One explanation for this could be 
that our method of sub-sampling modes in the wavevec- 
tor shells is not actually selecting representative samples 
of all the four-point configurations in each shell, which 
would be fixed by a more precise sub-sampling procedure. 

Because pursuing this algorithm further would signifi- 
cantly increase the computational complexity while lim- 
iting the range of applications, we instead turned to using 
multiple small-modes simulations with large-scale mode 



resampling to further explore the covariance matrix esti- 
mates as described in the next Section. 

3.2.3. Covariance matrix estimates from mode-resampling 

The power spectrum variance and correlation coeffi- 
cients estimated from 500 resamplings of the large-scale 
modes in a single small-mo des simulation are shown by 
the blue squares in Figure [lOl Note that we use the 



sa me binning in wavenu mber in all the cases shown 



as 



Takahashi et al. ( 2009 1 in order to make an accurate 
comparison. While the small-scale variance estimate is 
boosted beyond the Gaussian prediction, it is system- 
atically bias ed low relative to the results of |Takahashi| 
et al. ( 2009 ) . At intermediate wavenumbers just larger 
than fcthresh, the estimated variance is only ~ 20% of the 
true value. This is not too surprising because it is at 
these scales where our algorithm for adding large-scale 
modes is expected to be most inaccurate. This is be- 
cause our approximation that longer wavelength modes 
simply change the effective local matter density must be 
wrong when the longer-wavelengths are nearly equal to 
the smaller wavelengths to be perturbed. The modes 
just larger and smaller than fcthrosh may both contribute 
to coherent structures in the mass density distribution 
so that resampling only some of the Fourier modes in 
these structures would severely underestimate the true 
variance on these scales. Also, we should note that the 
modes with wavenumbers just smaller than /cthrcsh are 
likely not strictly Gaussian distributed in the all-modes 
simulation, so that their true distribution and coupling 
to small modes is not properly represented. 

The right panel in Figure 10 shows four rows of the 
matrix of estimated correlation coefficients of the power 
spectrum. It is significant that the mode resampling 
generates nonzero correlation coefficients between the 
large and small scales (blue crosses - lower panels es- 
pecially), which are identically zero in the unperturbed 
small-modes simulation. However, the correlation coef- 
ficients are again biased low by multiplicative factors of 
two or more. 

We have checked that the biases in the estimated power 
spectrum variance and correlation coefficients do not 
change when using more than 500 large-scale mode real- 
izations. Instead we speculate that this bias is real and 
indicates that (at least for fixed /cthrosh) the mode re- 
sampling algorithm generates realizations of the matter 
density from only a subset of the full nonlinear density 
probability distribution. This makes sense intuitively if 
we realize that the mode resampling as implemented here 
can never generate a new set of halos in the simulation 
box, but only moves and merges or de-merges the halos 
already present in the small-modes simulation. Several 
previous investigations have shown that the phase cor- 
relations in the Fourier modes of the density field in a 
dark matter simulation are dominate d by the masses and 



positions of the most mas sive halos (Watts et al. 2003 
Neyrinck & Szapudi 2007). It makes sense then that the 



power spectrum variance on small-scales would also be 
quite sensitive to the particular realization of the most 
massive halos. However, with a larger simulation volume 
the covariance estimates might be less dependent on the 
particular realization of the massive halos because many 
different halo configurations would be present. 
If our mode resampling algorithm is simply sampling 
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Figure 10. Left: Power spectrum variance normalized to that for a Gaussian density field with the same linear power. Right: Select 
correlation coefficients of the power spectrum covariance. Each panel corresponds to a row of the matrix of correlation coefficients, with 
the row labeled by the annotated k value in each panel (which is also where the correlation coefficients are all equal to one each panel). 
The vertical da shed lines show t he value of fcthrcsh = Skp. The colors and point styles are the same for both panels. The black circles 



{2009 ). The red diamonds are derived from the sample covariance estimate from 500 resamplings of each 
The blue squares are the estimates using 500 resamplings of only 1 small-modes simulation. Finally, the 



are values from [Takahashi et al, 

of 20 small-modes simulations. The blue squares are the estimates using 500 resamplings of only 1 small-modes simulation. Finally, 
green x's are the estimates when using 20 simulations without any resamplings (i.e., the sample covariance estimate that would be derived 
without knowledge of the algorithm in this paper.) 

a subset of the density probability distribution, then the 
bias in the power spectrum variance and correlation co- 
efficients should decrease if we use more than one simu- 
lation. To check this hypothesis we ran 19 more small- 
modes simulations (for 20 total) with different random 
number seeds in the initial conditions and added 500 re- 
samplings of the large-scale modes to each simulation 
(for a total of 10000 realizations of the density field). 
The sample covariance estimates impro ve c onsiderably 
as shown by the red diamonds in Figure [lOj The biases 
in both the variance and the correlation coefficients are 
much smaller over all scales. In addition, the extremely 
low dip in the variance estimates at intermediate scales 
when using only one simulation has disappeared entirely 
with 20 simulations. 

Finally, we show an estimate of the covariance matrix 
using 20 simulations without any resampling of large- 
scale modes with the green crosses in Figure 10 With 



We show the distribution in the variance estimates for 
wavenumber bins with centers in 0.295 < fc/(Mpc/fe) < 
0.395 with the boxes in the left panel of Figure [Tl] The 
horizontal lines in the centers of the boxes denote the 
median of the variance estimates and the box top and 
bottom denote the first and third quartiles. We have se- 
lected only high-fc modes to plot in Figure [Tl] because 
the variance estimate in this range has a roughly con- 
stant bias as seen in Figure [TO] and we are primarily in- 
terested in the accuracy of the covariance estimates on 
small scales where the deviation from the Gaussian model 
is most severe. We show how the normalization of the 
fit to the covariance dispersion changes with the number 
of mode resamplings per iV-body simulation in the right 
panel of Figure |11[ The reduction in the number of iV- 
body simulations needed to give a certain error in the 
estimate covariance flattens around 100 resamplings per 
simulation. 



such a small number of simulations the statistical noise 
dominates in the estimator. Comparing with the esti- 
mator with 500 resamplings of each of the A^-body sim- 
ulations shows how much our mode-addition algorithm 
reduces the statistical errors. 

3.2.4. Convergence rate of covariance estimates 

Although with 20 small-modes simulations the bias in 
the power spectrum covariance estimate is much smaller 
than with only o ne s imulation, a s mall bias is still no- 
ticeable in Figure [TO) In Figure [TT] we attempt to quan- 
tify the improvement in the covariance matrix estimates 
as more small-simulations are added by plotting the 
r.m.s. dispersion in elements of the power spectrum vari - 



Takahashi et al. 



ance relative to the result from |Takahashi et al. (2009) 



(20091 performed a similar analysis 
and found that the dispersion in their sample variance 
estimators was well fit by a — 2/Nr, where Nr is the 
number of iV-body simulations (or realizations) used to 
compute the estimator (shown as the dashed black line 
in Figure 11). We find a nonlinear least squares fit to 
our variance estimates of 2/(8. 57V °'^), which is shown 
by the solid red line in Figure |11[ From this we con- 
clude that the sample variance estimator using resam- 
pled large-scale modes can achieve similar accuracy as 
the standard sample variance estimator but with ~ 8 
times fewer TV-body simulations. However, it is diffi- 
cult to extrapolate this convergence rate far beyond the 
20 simulations we have actually performed because small 
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Figure 11. Left; Dispersion in components of the power spectrum covariance matrix estimates as functions of the numbe r of A'^-body 
simulations used to estimate the covariance. The dispersion is measured relative to the result from [Takahashi et ah] | |2009| using 5000 
simulations. Each A'^-body simulation has 500 resamplings of the large-scale modes, so a point on the plot has (number of A'^-body 
simualtions) X 500 total power spectra used to estimate the covariance. The boxes indicate the medians (central black lines) and quartiles 
(top and bottom of the boxes) for the dispersion of the diagonal covar iance elements within i: he wavenumber bins 0.295 < fc/(h/Mpc) < 
0.395. The dashed black line indicates the fit to the dispersion found in [Takahashi et al.| l |2009l l. The red line is a fit to the medians of our 
dispersion measurements. Right: Normalization of the fit to cr^ov versus Nr as a function of the number of mode resamplings per A'^-body 
simulation. That is, we plot the fit parameter A = 2/ (^cr^^^N^) just as for the red solid line in the left panel, but using between 1 and 500 
mode resamplings per A^-body simulation. (We fit a different b parameter for each value on the abcissa, but do not show these values.) The 
error bars denote the 95% confidence intervals on the nonlinear least-squares parameter estimate. The dashed grey line shows the expected 
scaling if each mode resampling was equivalent to running a new A^-body simulation. 



changes in our choice of wavenumber bins or in our least- 
squares fit can cause large extrapolation errors. 

4. DISCUSSION AND CONCLUSIONS 

We have demonstrated an algorithm to resample the 
large-scale Fourier modes of the density in an iV-body 
simulation that successfully reproduces the nonlinear 
power spectrum and significantly reduces the number of 
simulations required to estimate the power spectrum co- 
variance matrix. We expect our algorithm to aid the cal- 
culation of uncertainties for observational probes of the 
matter power spectrum such as weak lensing and galaxy 
clustering. Because we can quickly provide multiple re- 
alizations of the density field, it is possible to estimate 
uncertainties after applying appropriate survey windows 
to the density, which can have significant impacts on the 
estimated power spectrum covariances. This could be a 
distinct advantage over approaches that either precom- 
pute the power spectrum covar iance without consider- 
ing the survey window (f Sembolo ni et al.|2007 Takahashi 
et al.|200"9 Sato et aLp OOO ) , or use approximate or ana- 
lytic methods that may not hiclud e the survey windows 
correctly (e.g. Crocce et al. 2010). We have no reason 
to believe the extension to redshift space will pose any 
particular challenges. However, note that we have not 
yet tested the perturbation of the particle velocities so 
our results are limited to real-space statistics for now. 

Our algorithm uses both Zel'dovich displacements as 
well as time-perturbations to match the effective local 
linear growth in over- or overdense regions. It is intrigu- 
ing that neither of these operations alone is sufficient to 
reconstruct even the nonlinear power spectrum. Evaluat- 



ing the time-perturbation at the particle positions prior 
to the Zel'dovich move gives a bad result, indicating that 
the early-time movement of the particles is significant in 
determining the later growth of structures. Any exten- 
sions of our algorithm will likely have to consider both 
of these components carefully. 

Because applying the time-perturbation to remove 
large-scale modes from an existing simulation generally 
requires many simulation snapshots far into the future, 
we found it onerous (or even impossible if resampling 
too many modes in a box that is too small) to use ex- 
isting simulations for mode-resampling. However, if the 
simulation volume is at least several Gpc on a side, then 
according to Figure[7]it may be possible to subtract a lim- 
ited number of large-scale modes with a feasible amount 
of computation. For simulations volumes smaller than 
many Gpc on a side, we advocate running a new simula- 
tion with large-scale modes removed in the initial condi- 
tions, which can then have new large-scale modes added 
directly in later snapshots using our algorithm. Note 
that even for adding new large-scale modes the time- 
perturbation requires many snapshots closely spaced in 
scale factor both in the past and future of the desired 
time (but the simulation does not have to be run as 
far into the future as is needed for removing large-scale 
modes). This means that application of our algorithm 
for, e.g. parameter estimation from a galaxy survey, 
will require special-purpose A^-body simulations that will 
save significant computation time for the error analysis. 

Our algorithm has other applications such as includ- 
ing scaling of halo masses and concentrations so that 
mock galaxy catalogues could be constructed from each 
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of the resampled density fields. In our view, this is the 
only viable method to obtain robust uncertainties for es- 
timating cosmological parameters from galaxy clustering 
and BAOs. Our algorithm may even provide a feasible 



method to go beyond the work of Scfusatti et al. 



(20061 



for estimating the covariancc of higher-order correlations. 
For parameter estimation with weak lensing, our resam- 
pled density fields could be directly input to existing ray- 
tracing pipelines. However, one would need resampled 
densities for each of the lens planes (which could extend 
to z ~ 3 for upcoming survey requirements) . This would 
require saving even more closely spaced snapshots of the 
small-modes simulation(s). 

The sample covariance matrix estimated from a large 
number of resamplings of a single simulation has nonzero 
correlations between large and small scales that are in- 
troduced by our mode addition algorithm. But, the co- 
variance elements are biased low relative to the bench- 

which used 



mark result from Takahashi et al 
5000 TV-body simulations 



( |2009[ ), 

A possible explanation for 



this bias is that our method of resampling large-scale 
modes samples only a subset of the full probability dis- 
tribution for the nonlinear density field. We have shown 
that the bias decreases monotonically as a function of 
the number of TV-body simulations used to estimate the 
covariance, which is consistent with this explanation. If 
our mode-resampling algorithm were perfect at adding 
large-scale fluctuations, the bias in the covariance that 
we find would indicate that the distribution of the small- 
scale nonlinear density field is not entirely determined by 
the coupling with large-scale modes. Rather, the phases 
of the small-scale Fourier modes in the initial conditions 
would be important as well, which contradicts the con- 
clusion of 'Little et al. (1991) who found little impact on 
the final nonlinear density when the small-scale phases 
were randomized in the initial conditions. This could be 
because the eye notices mostly the large-scale structures 
in a dotplot of the particle positions but there can be sig- 
nificant small-scale deviations that are less noticeable (as 
shown by the "Zel'dovich only" lines in Figure 8 versus 
the "Zel'dovich only" lower left panel of Figure 6 1. 

One area for future investigation is to study the de- 
pendence of the bias in the power spectrum sample co- 
variance (or other measures of the distribution of the 
nonlinear density) with fcthrcsh- This would quantify the 
relative importance of different large-scale modes on the 
nonlinear density distribution. This could be a promising 
way to better understand the "loss" of information about 
cosmol ogy in the matter power spectrum on transhnear 



scales dRimes fc Hamiltmi 2006 Neyrinck et al. 2006 



iNeyrinck & Szapudi 2007p . One could also quantity oil 



which scales (it any] the density in a large simulation be- 
comes ergodic (so that the distribution of the density in 
sub- volumes of the simulation is equivalent to the distri- 
bution of the density field in an ensemble of large simula- 
tions) . The fact that our small-scale covariance estimates 
have such a large bias for /cthrosh ~ 0.5/i/Mpc (and in a 
single simulation) suggests that assuming ergodicity on 
these scales is a fallacy when defining estimators of, at 
least some, summary statistics of the density. 

Another application of our algorithm could be for tiling 
smaller simulation volumes to obtain the (properly corre- 
lated) density field in volumes t hat are too large to sim- 
ulate, as originally proposed by Cole (1997). This would 



probably be most effective if we run several of our "small- 
modes" simulations and then tile them together by re- 
sampling the modes both larger and smaller than the fun- 
damental frequency of the simulation volume. A related 
application is re-running resimulation studies (that take 
a sub-volume of high-resolution A^-body simulation and 
add additional physics) with different large-scale mode 
realizations. This could be a relatively fast way to ob- 
tain different realizations of the resimulations with only 
one high-resolution A^-body simulation at hand. 

Finally, we hope to apply our algorithm for obtain- 
ing estimates of both the mean power spectrum and 
its covariance as needed in simulation emulator frame- 



works ( [Heitmann et al.|2009[ [Schneider et al.||2011[ ). Be 

cause we can estimate the power spectrum covariance 
with an easily achievable number of simulations, it will 
now be feasible to build an emulator for the cosmol ogical- 
param eter dependence of the covariance (Schneider et al. 
2008|) . Note tha t this paper is complementary to Angulo 



& White (2010) who described how to rescale the out 
puts ot a single A^-body simulation to mimic the outputs 
with a different cosmology, but who did not consider the 
time-perturbations we use here. In combin ation, our al- 
gorithm and that oflAngulo & White ( 2010 1 may serve to 
drastically reduce the number of simulations needed to 
construct simulation emulators for cosmological analysis. 
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